home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / SQRT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-21  |  2.9 KB  |  177 lines
  1. /*                            sqrt.c
  2.  *
  3.  *    Square root
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, sqrt();
  10.  *
  11.  * y = sqrt( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the square root of x.
  18.  *
  19.  * Range reduction involves isolating the power of two of the
  20.  * argument and using a polynomial approximation to obtain
  21.  * a rough value for the square root.  Then Heron's iteration
  22.  * is used three times to converge to an accurate value.
  23.  *
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *
  29.  *                      Relative error:
  30.  * arithmetic   domain     # trials      peak         rms
  31.  *    DEC       0, 30        2000       2.1e-17     5.2e-18
  32.  *    IEEE      0,1.7e308   30000       1.7e-16     6.3e-17
  33.  *
  34.  *
  35.  * ERROR MESSAGES:
  36.  *
  37.  *   message         condition      value returned
  38.  * sqrt domain        x < 0            0.0
  39.  *
  40.  */
  41.  
  42. /*
  43. Cephes Math Library Release 2.1:  December, 1988
  44. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  45. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  46. */
  47.  
  48.  
  49. #include "mconf.h"
  50.  
  51. extern double SQRT2;  /*  SQRT2 = 1.41421356237309504880 */
  52.  
  53. double sqrt(x)
  54. double x;
  55. {
  56. int e;
  57. short *q;
  58. double z, w;
  59. double frexp(), ldexp();
  60.  
  61. if( x <= 0.0 )
  62.     {
  63.     if( x < 0.0 )
  64.         mtherr( "sqrt", DOMAIN );
  65.     return( 0.0 );
  66.     }
  67. w = x;
  68. /* separate exponent and significand */
  69. #ifdef UNK
  70. z = frexp( x, &e );
  71. #endif
  72. #ifdef DEC
  73. q = (short *)&x;
  74. e = ((*q >> 7) & 0377) - 0200;
  75. *q &= 0177;
  76. *q |= 040000;
  77. z = x;
  78. #endif
  79.  
  80. /* Note, frexp and ldexp are used in order to
  81.  * handle denormal numbers properly.
  82.  */
  83. #ifdef IBMPC
  84. z = frexp( x, &e );
  85. q = (short *)&x;
  86. q += 3;
  87. /*
  88. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  89. *q &= 0x000f;
  90. *q |= 0x3fe0;
  91. z = x;
  92. */
  93. #endif
  94. #ifdef MIEEE
  95. z = frexp( x, &e );
  96. q = (short *)&x;
  97. /*
  98. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  99. *q &= 0x000f;
  100. *q |= 0x3fe0;
  101. z = x;
  102. */
  103. #endif
  104.  
  105. /* approximate square root of number between 0.5 and 1
  106.  * relative error of approximation = 7.47e-3
  107.  */
  108. x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
  109.  
  110. /* adjust for odd powers of 2 */
  111. if( (e & 1) != 0 )
  112.     x *= SQRT2;
  113.  
  114. /* re-insert exponent */
  115. #ifdef UNK
  116. x = ldexp( x, (e >> 1) );
  117. #endif
  118. #ifdef DEC
  119. *q += ((e >> 1) & 0377) << 7;
  120. *q &= 077777;
  121. #endif
  122. #ifdef IBMPC
  123. x = ldexp( x, (e >> 1) );
  124. /*
  125. *q += ((e >>1) & 0x7ff) << 4;
  126. *q &= 077777;
  127. */
  128. #endif
  129. #ifdef MIEEE
  130. x = ldexp( x, (e >> 1) );
  131. /*
  132. *q += ((e >>1) & 0x7ff) << 4;
  133. *q &= 077777;
  134. */
  135. #endif
  136.  
  137. /* Newton iterations: */
  138. #ifdef UNK
  139. x += w/x;
  140. x = ldexp( x, -1 );    /* divide by 2 */
  141. x += w/x;
  142. x = ldexp( x, -1 );
  143. x += w/x;
  144. x = ldexp( x, -1 );
  145. #endif
  146.  
  147. /* Note, assume the square root cannot be denormal,
  148.  * so it is safe to use integer exponent operations here.
  149.  */
  150. #ifdef DEC
  151. x += w/x;
  152. *q -= 0200;
  153. x += w/x;
  154. *q -= 0200;
  155. x += w/x;
  156. *q -= 0200;
  157. #endif
  158. #ifdef IBMPC
  159. x += w/x;
  160. *q -= 0x10;
  161. x += w/x;
  162. *q -= 0x10;
  163. x += w/x;
  164. *q -= 0x10;
  165. #endif
  166. #ifdef MIEEE
  167. x += w/x;
  168. *q -= 0x10;
  169. x += w/x;
  170. *q -= 0x10;
  171. x += w/x;
  172. *q -= 0x10;
  173. #endif
  174.  
  175. return(x);
  176. }
  177.